/*!  \file  SteMR.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      06/29/2013 11:01:22 PM
 *   \copyright GNU Public License.
 */
#include "SteMR.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void dstemr_(
const char   * jobz, 	
const char   * range, 	
const int    * n,
      double * d,
      double * e,
const double * vl,
const double * vu,
const int    *	il,
const int    *	iu,
      int    *	m,
      double * w,
      double * z,
const int    *	ldz,
const int    *	nzc,
      int    *	isuppz,
      int    *	tryrac,
      double * work,
const int    * lwork,
      int    * iwork,
const int    * liwork,
      int    * info);
}
SteMR::SteMR(): TriEigen("Lapack dstemr") {	
	setCompileTime(__DATE__, __TIME__); 
	canBeParallel = false;
}
bool SteMR::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
void SteMR::computeRange( const Matrix & T, Matrix & Sigma, Matrix & Z, 
	const int indexLow, const int indexUp) {
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Z.isOrthogonal());
	char jobz = 'V', range = 'I';
	int n = T.cols(), info = -1, ldz = Z.lda();
	double vl = 0, vu = 1; int il = indexLow, iu = indexUp;
	assert(il >= 1 && iu >= 1 && il <= n && iu <= n);
	int m = 0; // total number of eidgen values found
	double * w = NULL; allocate(w, n);//Sigma.getModData();
	double * z = Z.getModData(); 
	if (Z.isRowMajor()) Z.makeColMajor();
	double * d = NULL, *e = NULL, *work = NULL; allocate(d, n); allocate(e, n);
	memcpy(d, T.getData()       ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	e[n-1] = 0;	
	int tryrac = 1;
	int nzc = iu - il + 1;
	int * isupp = NULL, *iwork = NULL;
	allocate(isupp, 2 * nzc);
	double optW = 0; int optIW = 0, lwork = -1, liwork = -1;
	dstemr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &m, w, z, &ldz,
			  &nzc, isupp, &tryrac, &optW, &lwork, &optIW, &liwork, &info);
	lwork = (int) optW; liwork = optIW;
	lwork  = max( lwork, max(jobz == 'V' ? 18 * n : 12 * n, 1));
	liwork = max(liwork, max(jobz == 'V' ? 10 * n :  8 * n, 1));
	allocate(work, lwork); allocate(iwork, liwork);
	if (work && iwork) {
		dstemr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &m, w, z, &ldz,
		  &nzc, isupp, &tryrac, work, &lwork, iwork, &liwork, &info);
		success &= info == 0;
	}
	delocate (iwork);	delocate ( work);	delocate (isupp);
	delocate ( e   );	delocate ( d   ); delocate (w);
}
void SteMR::compute(const Matrix & T, Matrix & Sigma, Matrix & Z ){
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Z.isOrthogonal());
	char jobz = 'V', range = 'A';
	int n = T.cols(), info = -1, ldz = Z.lda();
	double vl = 0, vu = 1; int il = 0, iu = n-1;
	int m = 0; // total number of eidgen values found
	double * w = Sigma.getModData();
	double * z = Z.getModData(); Z.makeColMajor();
	double * d = NULL, * e = NULL, *work = NULL;
	allocate(d, n); allocate(e, n);
	memcpy(d, T.getData()    ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	e[n-1] = 0;
	int nzc = max(n,1);
	int tryrac = 1;
	int * isupp = NULL, *iwork = NULL;
	allocate(isupp, 2 * n);
	double optW = 0; int optIW = 0, lwork = -1, liwork = -1;
	dstemr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &m, w, z, &ldz,
			  &nzc, isupp, &tryrac, &optW, &lwork, &optIW, &liwork, &info);
	lwork = (int) optW; liwork = optIW;
	lwork  = max( lwork, max(jobz == 'V' ? 18 * n : 12 * n, 1));
	liwork = max(liwork, max(jobz == 'V' ? 10 * n :  8 * n, 1));
	allocate(work, lwork); allocate(iwork, liwork);
	if (work && iwork) {
		dstemr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &m, w, z, &ldz,
		  &nzc, isupp, &tryrac, work, &lwork, iwork, &liwork, &info);
		success &= info == 0;
	}
	delocate (iwork);	delocate ( work);	delocate (isupp);
	delocate ( e   );	delocate ( d   );
};
BiNESteMR::BiNESteMR() : BiSVD( "NE dstemr" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteMR()));
	setCompileTime(__DATE__, __TIME__); 
};
BiNE2SteMR::BiNE2SteMR() : BiSVD( "NE2 dstemr" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteMR()));
	setCompileTime(__DATE__, __TIME__);
};
BiGKSteMR::BiGKSteMR() : BiSVD( "GK dstemr" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteMR()));
	setCompileTime(__DATE__, __TIME__);
};
REGISTER(SteMR);
REGISTER(BiNESteMR);
REGISTER(BiNE2SteMR);
REGISTER(BiGKSteMR);
